###### Script - Ecosystems and Ordering: a dataset  
### transforming "dataset_ecosystem_wide.xlsx/csv" in "dataset_ecosystem_long.xlsx/csv"

#set directory
setwd()

#libraries
library(tidyverse)
library(readxl)
library(writexl)
library(sf)
library(maps)

dataset <- read_excel("dataset_ecosystem_wide.xlsx")
acronyms <- read_excel("Cooperation Initiatives.xlsx")

d <- dataset

dd<-c()
for(r in 1:nrow(d)){
  dd<-rbind(dd,d[r, c(1:11,12:22)])
  
  ddd<-d[r, c(1:11,23:33)]
  colnames(ddd)<-colnames(dd)
  dd<-rbind(dd,ddd)
  
  ddd<-d[r, c(1:11,34:44)]
  colnames(ddd)<-colnames(dd)
  dd<-rbind(dd,ddd)
  
  ddd<-d[r, c(1:11,45:55)]
  colnames(ddd)<-colnames(dd)
  dd<-rbind(dd,ddd)
  
  ddd<-d[r, c(1:11,56:66)]
  colnames(ddd)<-colnames(dd)
  dd<-rbind(dd,ddd)
  
}


dataset_ecosystem_long <- dd %>% 
  unique() %>% 
  filter(ecosystem_number_adjancent_countries >= 4) %>% 
  filter(!(ecosystem_regional_cooperation %in% c(1,2) & is.na(cooperation_1)))

# Visualization: Figure 4
acronyms <- read_excel("ecosystem_cooperation_initiatives_acronyms.xlsx")

dataset_ecosystem_long %>%
  left_join(acronyms, by = c("cooperation_1" = "cooperative_initiatives2")) %>%
  filter(!is.na(acronyms)) %>% 
  group_by(cooperation_1, cooperation_1_type, acronyms) %>%
  summarize(n = n(), .groups = 'drop') %>%
  ggplot(aes(y = fct_rev(acronyms), x = as.factor(cooperation_1_type))) +
  geom_point(aes(size = n)) +
  labs(x = "Cooperation Type", y = "", size = "Frequency") 

#long file
write_xlsx(dataset_ecosystem_long, "dataset_ecosystem_long.xlsx")
write_csv(dataset_ecosystem_long, "dataset_ecosystem_long.csv")


####### maps
#shapefiles
terrestrial_shapefile <- st_read("official/wwf_terr_ecos.shp") ## available at: https://www.worldwildlife.org/publications/terrestrial-ecoregions-of-the-world
marine_shapefile <- st_read("MEOW/meow_ecos.shp") ## available at: https://www.worldwildlife.org/publications/marine-ecoregions-of-the-world-a-bioregionalization-of-coastal-and-shelf-areas
freshwater_shapefile <- st_read("GIS_hs_snapped/feow_hydrosheds.shp") ## available at: https://www.feow.org/download
world_borders <- map('world', fill = TRUE, col = 'white', plot = FALSE)


# Visualization: Figure 1
terrestrial_dataset <- dataset %>% 
  filter(teows_meows_feows_provinces == "teow")
terrestrial_dataset$ecosystem <- trimws(terrestrial_dataset$ecosystem)
terrestrial_dataset$ecosystem <- tolower(terrestrial_dataset$ecosystem)

terrestrial_shapefile$ECO_NAME <- trimws(terrestrial_shapefile$ECO_NAME)
terrestrial_shapefile$ECO_NAME <- tolower(terrestrial_shapefile$ECO_NAME)

merge_terrestrial <- merge(terrestrial_shapefile, terrestrial_dataset, by.x = "ECO_NAME", by.y = "ecosystem")
merge_terrestrial$ecosystem_regional_cooperation <- as.factor(merge_terrestrial$ecosystem_regional_cooperation)

ggplot() +
  geom_polygon(data = world_borders, aes(x = long, y = lat, group = group), color = "gray50", 
               fill = NA, alpha = 0.5) +
  geom_sf(data = merge_terrestrial, aes(fill = ecosystem_regional_cooperation), alpha = 0.7) + 
  scale_fill_manual(values = c("lightgoldenrod1", "green3", "darkgreen"),
                    na.value = "white") +
  theme_minimal() +
  theme(legend.title = element_blank(), 
        legend.text = element_text(size = 12),
        axis.text = element_blank(),
        axis.title = element_blank()) 


# Visualization: Figure 2
marine_dataset <- dataset %>% 
  filter(teows_meows_feows_provinces == "meow")
marine_dataset$ecosystem <- trimws(marine_dataset$ecosystem)
marine_dataset$ecosystem <- tolower(marine_dataset$ecosystem)

marine_shapefile$ECOREGION <- trimws(marine_shapefile$ECOREGION)
marine_shapefile$ECOREGION <- tolower(marine_shapefile$ECOREGION)

merge_marine <- merge(marine_shapefile, marine_dataset, by.x = "ECOREGION", by.y = "ecosystem")
merge_marine$ecosystem_regional_cooperation <- as.factor(merge_marine$ecosystem_regional_cooperation)

ggplot() +
  geom_polygon(data = world_borders, aes(x = long, y = lat, group = group), color = "gray50", 
               fill = NA, alpha = 0.5) +
  geom_sf(data = merge_marine, aes(fill = ecosystem_regional_cooperation), alpha = 0.7) + 
  scale_fill_manual(values = c("lightgoldenrod1", "darkorchid3", "royalblue1"),
                    na.value = "white") +
  theme_minimal() +
  theme(legend.title = element_blank(),
        legend.text = element_text(size = 12),
        axis.text = element_blank(),
        axis.title = element_blank()) 

# Visualization: Figure 3
freshwater_dataset <- dataset %>% 
  filter(teows_meows_feows_provinces == "feow")
freshwater_dataset$ecosystem_id <- trimws(freshwater_dataset$ecosystem_id)
freshwater_dataset$ecosystem_id <- tolower(freshwater_dataset$ecosystem_id)

freshwater_shapefile$FEOW_ID <- trimws(freshwater_shapefile$FEOW_ID)
freshwater_shapefile$FEOW_ID <- tolower(freshwater_shapefile$FEOW_ID)

merge_freshwater <- merge(freshwater_shapefile, freshwater_dataset, by.x = "FEOW_ID", by.y = "ecosystem_id")
merge_freshwater$ecosystem_regional_cooperation <- as.factor(merge_freshwater$ecosystem_regional_cooperation)

ggplot() +
  geom_polygon(data = world_borders, aes(x = long, y = lat, group = group), color = "gray50", 
               fill = NA, alpha = 0.5) +
  geom_sf(data = merge_freshwater, aes(fill = ecosystem_regional_cooperation), alpha = 0.7) + 
  scale_fill_manual(values = c("lightgoldenrod1", "darkorange1", "firebrick4"),
                    na.value = "white") +
  theme_minimal() +
  theme(legend.title = element_blank(), 
        legend.text = element_text(size = 12),
        axis.text = element_blank(),
        axis.title = element_blank()) 
